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STABLE IMPLICIT AND EXPLICIT NUMERICAL METHODS FOR 


INTEGRATING QUASI -LINEAR DIFFERENTIAL EQUATIONS 
WITH PARASITIC -STIFF AND PARASITIC -SADDLE 
EIGENVALUES 
By Harvard Lomax 
Ames Research Center 

SUMMARY 


Certain classes of coupled, quasi -linear, ordinary, differential equa- 
tions contain eigenvalues in their associated matrix which make them difficult 
to integrate hy means of conventional numerical differencing schemes, even 
when the solutions are continuous and nonsingular. Two classes of such 
"parasitic" eigenvalues are defined and general ways in which their effects 
can he suppressed are discussed. 


INTRODUCTION 


Special classes of differential equations ("both ordinary and partial) , 
termed quasi-linear, are defined by the condition that the highest order deriv- 
ative terms appear explicitly and to the first power only. For ordinary dif- 
ferential equations these can be written 

[B]w' = C (l) 

where the elements of [B] and C can depend upon w and upon the independent 
variable t. This equation can be written formally 

w» = [B]" 1 C = F(w,t) ' (2) 

Cases in which [B] -1 does not exist (det(B) = 0) are of special interest and 
lead to the study of critical points, and in particular, for the discussion in 
this paper, to the study of saddle points. 

— > 

We are principally concerned -with the practical situation when F has a 
nonlinear dependence upon w- However, we assume this dependence is continu- 
ous at least through the first derivative- In such cases, we can make a local 
Taylor series expansion of F about some reference point n, where t = nh 
and h is a small (step) interval. Equation (2) then becomes 



f 


where [A n ] is the Jacobian of F with respect to w and 

f n = ^n " [A n ]w n (3h) 

and er-(- is bounded as h -*■ 0. 

If the term h 2 er^ is neglected in the interval nh ^ t ^ h(n + l) the 
result is, in that step, a set of coupled; ordinary; linear differential equa- 
tions. Further; if equations (2) are autonomous (i.e.; F does not depend 
explicitly on t); they are linear equations with constant coefficients. In 
carrying out practical numerical calculations; the elements of F (or [An] and 
%> if that form is used) are varied from step to step. This "captures" the 
nonlinear effects by embedding local polynomials of O(h^) in the calculations 
as they proceed. If equations (3) are used; w is approximated to 0(h 3 ). 
Further; if in each successive step a method is chosen that is stable for the 
local linearized form in that step; overall or "global" stability is assured 
in the sense described below. 

There are two common approaches to the numerical integration of differen- 
tial equations . One is to use explicit differencing formulas and apply them 
directly to equations (2); without ever actually formulating equations (3)* 
These methods can have unlimited accuracy (since the local linearization is 
never carried out) and; in the author's experience; have stability properties 
that correlate extremely well with the eigenvalues of [A n ] (as if the local 
linearization actually had been performed). However; all explicit methods 
have a finite (and; practically speaking; rather limited) stability boundary. 
The optimization of explicit methods; from the point of view of stability; is 
treated in reference 1; and briefly reviewed in the section Highly Stable 
Explicit Methods. 

Another approach is to use implicit methods. A general discussion of 
this approach is the intended contribution of this paper. In the implicit 
case, equations (2) are in fact put in the form of equations (3) and treated 
as coupled linear equations in a single step. The elements of [A n ] must be 
reevaluated at each step, although it is consistent with the accuracy involved 
to do this numerically rather than analytically. For example, use of the 
formula 

_ dFi F i (l.01 Wj) - F ± (0.99 w-j ) 
a ij = dvn * ~ " 0.02 w j 

has given good practical results. Since calculating all the a±j for each 
step can be quite time-consuming, it is advisable to know when it is and when 
it is not worthwhile. This leads directly to the subject discussed next. 


GENERAL COMMENTS ON THE CONSTRUCTION OF NUMERICAL METHODS 
FOR INTEGRATING ORDINARY DIFFERENTIAL EQUATIONS 


It is unlikely that "new" combinations of linear equations that connect a 
function, u, and its derivative, u', at a series of reference points, 
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equispaced or not , will improve existing methods for the numerical integration 
of general sets of coupled ordinary differential equations having the form 
given hy equations (2). If nothing special is known a priori about the dif- 
ferential equations, the standard fourth order Runge-Kutta method is probably 
the "best." It is self -starting, has low storage capabilities, is easy to 
program, has good accuracy 0(h 5 ), and, as we shall presently see, is more 
stable than any of the standard predictor -corrector processes (e-g. , Hamming f s, 
Adams -Moult on , etc.). 

Nevertheless, it is still popular to publish numerical integration meth- 
ods for general classes of ordinary differential equations (see, e.g., 
refs. 2-5). These methods are meant to compete with, or improve upon, the 
classical ones mentioned above. Under what qualifications they do, or do not, 
is a subject discussed in reference 6. For the most part, their development 
seems to be motivated by a popular misconception that one can "break through" 
the Dahlquist stability "barrier." Now the Dahlquist stability theorem is the 
result of a rigorous proof (see ref. 7) based upon certain premises. One of 
these premises is that the function and its derivative are evaluated at, and 
only at, a series of equispaced points. Another is that a predictor is fol- 
lowed by only one corrector. Under these conditions Dahlquist shows that no 
linear stable equations relate u n +q and u^+q, i = 1, 2, • • •, k (i.e., k-step 
methods) with an accuracy of order greater than k + 2 , if k is even, or 
k + 1, if k is odd. But if the conditions described above are violated in 
any way, the theorem no longer applies. 

In order to avoid the confusion and ambiguity that can arise when we try 
to classify explicit methods according to (computational) step size and step 
number, the representative step size is introduced such that: 

H = the distance a solution is advanced after 

two evaluations of the derivative. (4) 


If both the accuracy and stability of a method are referenced to this parame- 
ter, truly significant comparisons can be made amongst all types of numerical 
integration procedures involving linear connections between the function and 
its derivative. Let us next consider two examples. 


The standard, fourth order, Runge-Kutta method is usually referred to as 
a one -step method. When it is cast in predictor-corrector terminology, it can 
be written 


n+ F 

( 2 ) 

U i 


= u- 


u 


(3) 


'n+i 


u 


n+x 


■n 


= u, 
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- u- 
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= u n + l ~ 


1 

u n 
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+ 2u 
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+ U- 
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(5b) 

(5c) 

(5d) 
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Graphically (see definition of symbols in table i) it appears as shown in 
sketch (a) . From the sketch we see that the reference step size, H, is half 
that usually used in the computations. The example shows clearly that 


TABLE I. - DEFIHITIOUS OF SYMBOLS USED IE SKETCHES 


Symbol 

Represents 


□ 

Value of function predicted from data 
previous steps. 

at 

+ 

Derivative calculated using □ . 


o 

1. Value of corrected function using 
previous steps and predicted data at 
step, or 

2 . Value of function calculated using 
method . 

data at 
this 

implicit 

X 

Derivative calculated using O • 




Eq. (5a) 



Eq (5b) 


references to accuracy and stability on 
the basis of a calculation step size, 
h, and step number, k, (for eqs. (5) is 
k equal to 2 or 1?) are untrustworthy 
since the step designation is arbitrary. 
On the basis of H, equations (5) have 
the error terms 



Sketch (a).- Standard, fourth-order, Runge-Kutta 

method . for complementary and particular 

solutions (respectively) to the representative equation 


u y 


CJu + e 


JJLt 


( 7 ) 


There is, generally, no connection between step number, k, and the accu- 
racy and stability of a method. This has been discovered by various authors 
(see refs. 2-5), and has led to the introduction of terms such as "hybrid" and 
"combined" methods, terms which indicate (simply) that the premises required 
for the validity of Dahlquist's theorem have been abandoned. 

Typical of the hybrid or combined methods (for a complete analysis see 
ref. 6, pp. 93-95) is 
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(8a) 


u 


( 1 ) 


n+i .8 


u 


'n+p 


u + 
n+i 


= u , + 

n+i 


% ( 2S8u ' 


+ 22u’ - 230u^ 
n+i n n+o 


lV ) 

+o. ay 


, (i)' 

u + 25u v / -25u , 

n n+i.s n+o 


1)1 ) 

+o.sy 


(8b) 


which has the graphical representation shown in sketch (b) . Its error, when 

applied to equation (7), is given by 



<i 


-Wi 


4 + 


_L 




er A 


er 


720 

1 


(CTH) 5 

(m-H) 5 


(9a) 


(9b) 


"Predictor" 
Eq. (8a ) 


"Corrector" 
Eq (8b) 


'M- 720 

Clearly, equations (9) are much more 
accurate (when aH is small enough 
for the first missed term in a Taylor 
series expansion to give a valid 
approximation of the total error) than 
equations (5); but they pay the usual price: they are less stable. 


Sketch (b).- A typical "hybrid” or "combined" 
method . 


Now the stability of a method depends upon the roots to its characteris- 
tic equation when that method is used to difference equation ( 7 )- The stabil- 
ity of Runge-Kutta methods is analyzed in many places (see ; e.g., ref. 6 ). 

The characteristic polynomial for equations (5) (in terms of the displacement 
operator E = e ^(d/dx)) i s 


P(E) = E - l - oh - | (ah) 2 - \ (ah) 3 - ^ (ah) 4 ( 10 ) 

It has only one root (the principal one) and this root represents e a ^ 
exactly through the fourth order. The P(E) for the ith order Runge-Kutta 
method represents e ^* 1 exactly through the ith order , at which term it ends. 
The detailed behavior of these roots , as ah ranges from zero through a vari- 
ety of complex values, is shown on page 83 of reference 6 for the second, 
third, and fourth order methods. The stability boundary | crH | 0 (defined pre- 
cisely in the next section) is shown in 
sketch (c). From this sketch we see 
that the fourth order Runge-Kutta 
method is the best of all Runge-Kutta 
methods (i.e., minimizes machine com- 
puting time) when applied to coupled 
ordinary differential equations with 
parasitic (stiff) eigenvalues (a term 
defined precisely in the next section) . 

The characteristic equation for 
the method defined by equations (8) is 
a cubic, so that it contains two spuri- 
ous roots. However, both of these roots go to zero when h goes to zero. 
Because of this the method is said to have "Adams -Moult on stability" (this 



Sketch (c).- General stability boundary of 
Runge-Kutta methods. 
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property is characteristic of all Mams -Moulton methods), and equations (8) 
have guaranteed stability for small enough h. The detailed behavior of all 
the roots is shown on page 96 of reference 6, and the general stability bound- 
ary is shown in sketch (d). Although the method is quite accurate, it is not 

suitable for coupled equations with 
parasitic eigenvalues because of its 
low stability boundary. 

An analysis of the above examples 
(and many others presented in ref. 6) 
leads to the conclusion expressed in 
the first sentence of this part. This 
does, not mean, however, that valuable 
numerical methods (of the type being 
considered) cannot still be formulated. What appears to be needed are studies 
of special methods designed for special classes of equations (among which the 
quasi-linear equations ( 2 ) are already a special case). For example, if some- 
thing is known about the eigenvalues in the local associated matrix in equa- 
tions (3)> particular numerical methods can be constructed which are superior 
to the classical ones for the particular parameters involved. Such cases are 
considered below. 



Sketch (d).- General stability boundary for 
equations (8). 


DEFINITIONS AND TERMINOLOGY 


We wish now to introduce some general notation and define precisely some 
terms used in the later discussion. First, consider M coupled, linear, 
ordinary differential equations with constant coefficients and (for simplicity, 
though it is not essential) distinct eigenvalues. Then 


w- « £ = [A]w + f 
dt 


(id 


where, if c m j are constants dependent on the initial conditions, t = nh, and 
are the eigenvalues of [A], 

M 


w 


nm 


= ^ c mj( e<J ^) n + Articular S. , 


m = 1, 2, 


, M (12) 


J=i 


Next choose any set of linear, difference -differential equations with constant 
coefficients (e.g., Runge-Kutta, Adams -Moulton, Hamming’s, etc.) and symbolize 
it by the operator L. By itself L represents any linear set of differenc- 
ing operations composing a method; on the other hand, L e refers specifically 
to linear explicit methods, and Lj. to linear implicit ones. Thus if L e 


represents the predictor 


(1) . , 

u ' , ' - u + hu* 
n+i n n 


followed by the corrector 


u = u + 
n+i n 


2 h ( u »( + U a) , the operation L e (u f = cru + f) results in the 


difference equations 



u 


(1) _ 


n+i 


= (l + crh)u n + hf 


L n 


u 


n+i 


u n 4- J ah ( + u 


n+i 


n^ + 2 ^(^n+i + ^n) 


having the characteristic polynomial 


P(E) = E- l-crh-| cr 2 h 2 = 0 


For arbitrary L the operation L(u f = cru + f) results in a set of difference 
equations , the solution to which can always he written 

(13) 

where are the roots to the characteristic polynomial P(E) =0. The 

value of these roots depends on the choice of L and, in particular, 


u n 


Bi(Ai) n 


+ P.S. 


1=1 


Ai = gi(crH) , i = 1, 2, • • *, k 


(14a) 


Ai = e aE + 0(H Z ) (l4b) 

7\j_ are the spurious roots of the numerical method if i > 1 
are constants dependent upon the initial conditions 
H is the representative step size defined in equation (4) 

It should he noted that: 


The value of l depends upon the order of 
the local Taylor series which L embeds in 
the calculations. In general, l is inde- 
pendent of k even for stable L. Dahlquist 
imposes special conditions under which l 
and k are connected if stability is imposed. 


The function gj_ depends entirely upon the 
choice of L. 


(16) 


Now it can be shown (see ref. 6) that the operation L(w = 
in a set of coupled difference equations, the solutions to 
written 

M k 

^ Cmj(^ji) n + P -S- , a = 1; 2, 
j=i i=x 



[A]-w + f) results 
which can always he 


• • • , M (17) 
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"where 

Aji = gl (tfjH) , i = 1, 2, • • -,k ; j = 1, 2, • • •, M ( 18 a) 
Aj! = e°J E + 0(H Z ) (l8b) 


A.. are the spurious roots of the numerical method if i > 1 

Cm - are constants dependent upon the initial conditions 
J 

It should he noted that: 


If, and only if, 1 the same L is applied 
in a given calculation step, h, to all the 
equations in a coupled set , the numerical 
accuracy and stability do not depend upon 
the elements in [A n ] except as those 
elements affect the eigenvalues cf j . 

Under the same conditions as in (l9)j the 
functions g^ in equations (l4a) and (l8a) 
are identical. 

The following definitions can now be formulated. 

A set of differential equations is inherently stable if 
Re(tfj) s: 0, j = 1, 2, • • M. 


(19) 


( 20 ) 


( 21 ) 


If a set of differential equations is inherently stable,^ 
the set of difference equations formed by L(w =■ [A]w + f) 
has an induced instability if |Aji| = |gi(°jH)| > l, 

J = 1, 2, • • •, M; i = 1, 2, • • -,k. (Note that i = 1 
is included . ) 

If Re(crjH) £ 0, the value of crjH for which any increase 
in H makes |Aji| = |gi(ojH) | > 1 is labeled |ojH| c and 
called the general stability boundary, or simply, the 
stability boundary. 


( 22 ) 


(23) 


If we set Oj = tfje 10 (where <Jj and 9 are real), and let 
CTjH <0, the value of 0jH for which any increase in H 
makes |Ajp| = | gi (ffH) | > 1 is labeled |o&| c and called 
the real stability boundary. Similarly, the value of 
itf^H for which any increase in H makes 
|A ji | = |gp(iffjH) | > 1 is labeled |ibH| c and called the 
imaginary stability boundary. 


(24) 


1 Proof of ( 19 ) is given in reference 6. It should be noted that roundoff 
effects (which may be important) are neglected in all statements made in this 
report . 
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Methods for which | crH | c < 00 a ^ e called conditionally stable. 


( 25 ) 


All explicit methods are conditionally stable. A proof is 
given in reference 1. 

Certain classes (containing any order of accuracy) of 
implicit methods are unconditionally stable. This is 
discussed in the last section. 

The Oj can be divided into two classes: 

1. Driving eigenvalues,, q of them with subscript d, say, and 

2. Parasitic eigenvalues, M-q of them with subscript p 


(2 6 ) 

( 27 ) 


and the parasitic eigenvalues can be subdivided into two groups: parasitic- 

stiff and parasitic-saddle. The definitions follow: 


Driving combined with parasitic- 
stiff eigenvalues result when 


« C T 


M 


/ O' 

Cj( e J ) d » 


J=1 


0=1 



(28) 


Driving combined with parasitic- 
saddle eigenvalues result when 


Re(a) p > 0 
Re(a) d < 0 


I 


j=i 


/ O jH» 
i(e J ) c 


i 0 


( 29 ) 


M 

c j( eCJ ^ H )p = 0 (analytically) 

j=q. 
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THE CONTROL OF PARASITIC EIGENVALUES 


Consider the following numerical problem: How can one integrate coupled, 

quasi -linear, differential 'equations with both driving and parasitic eigen- 
values, letting the step size be determined by the driving ones, and, at the 
same time, introducing no significant errors due -to the parasitic ones? 

Of course, one could always locally linearize the equations and uncouple 
them. But this, in effect, requires not only finding the elements of [A], but 
also the eigenvalues and, what is more drastic, the eigenvectors themselves. 

The remarkable thing is that , although there are cases for which finding the 
elements of [A] is advisable, actually calculating the eigenvalues or eigen- 
vectors is unnecessary. It is a direct consequence of (19) that, under the 
qualifications noted, any standard, linear, differencing scheme integrates the 
coupled equations by seeking out each individual eigenvalue and integrating it 
as if the others did not exist (see footnote l) . Under these conditions, a 
step size can be chosen that will integrate one group of eigenvalues completely 
inaccurately but another group with high accuracy, and if the two groups are 
coupled together, the high accuracy would remain for the one in spite of the 
large errors in the others. (An example is given in ref. 8, pp. 35 - 36.) How- 
ever, in general, all the coupled solutions would be inaccurate, and the highly 
accurate results could only be recovered by uncoupling the final answers. 

While the last remark is interesting, it is not quite at the heart of the 
matter. Consider instead the following logic which is, 'in fact, the "solution" 
to the problem posed at the beginning of this section. 

1. Let the coefficients cj of (e^J^Jp in equations (28) and (29) be 
made small with respect to the coefficients of (e ff jH)^ by choosing appropriate 
initial conditions at the commencement of the numerical integration. 

2. Choose any numerical method, L, that is stable for all (ajH)p with 
total disregard as to its accuracy for them . 

3. Let the same L be both stable and accurate for all (djH)^. 

4. Then, in applying L to a set of differential equations in which 
(ffjH)p and (ffjH)a are coupled in any fashion, all of the coupled solutions 
will be accurately resolved and will represent the solution due to the driving 
eigenvalues as if the parasitic ones had been removed. 


HIGHLY STABLE EXPLICIT METHODS 


Some highly stable explicit methods for integrating coupled, autonomous , 
quasi-linear, differential equations with real, parasitic-stiff eigenvalues 
were developed in reference 1. For the sake of completeness, the subject is 
briefly reviewed in this section. 
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The characteristic equation for any explicit method; L e , is a monic 
polynomial of the general form 


+ P ((TH)^ -1 + • • • + P q (cJH) = 0 

A, " 1 


( 30 ) 


•where the Pj(ffH) are themselves polynomials in crH. Highly stable explicit 
methods for real; parasitic-stiff eigenvalues can he constructed; if L e is 
chosen so that the roots to its characteristic polynomial; equation ( 30 ); 

behave like those shown in sketch (e). 




A _ 
H°-h| c 


Sketch (e) . 


The principal root must correspond to 
the expansion of e 0 ^ through order 
Z if the accuracy of L e is to be 
O(H^). Aside from this; it is only 
necessary to make the magnitude of Ai 
and all the spurious roots Ai; i > 1 
(if there are any) less than one for 
as large a range of -crH as possible. 

In reference 1 one -root methods 
were developed for autonomous equa- 
tions by finding the b n in 


N 

Ai = 1 + tfh + | (ah ) 2 + ^ 'b n (oh) n = e ah + 0 (h 3 ) ( 3 l) 

n=3 


such that |cxh| c in sketch (e) was maximized in a least squares sense. (The 
actual optimum polynomials are unknown for N > 3-) A summary of the results 
is contained in the following chart. 


N 

k 

5 

Me 

12 

17 

Me 

5-5 

7-2 


6 

25 

8.6 


7 

8 

9 

10 

35 

45 

57 

70 

10.1 

11.4 

12.7 

14.1 


Notice that methods can be constructed that are highly stable for complex 
or even imaginary ay Furthermore; they can be constructed such that they 



Sketch (f). 


are stable in certain bands of 0 H 
and unstable in other bands. For exam- 
ple; taking again the case for real a 
(not only for simplicity; but since 
this is a practical situation) ; the 
b n in equation ( 3 l) could be chosen 
so that the curve in sketch (f) would 
result. The method L e which had 
such a (single) root in its 
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character! sti c polynomial would "be stable when used to integrate differential 
equations with rather large negative eigenvalues , provided the intermediate 
smaller set did not appear in the coupled group. This kind of method has not 

been thoroughly exploited., but it is 
not impractical. An example of its 
use is contained in the method sug- 
gested by Treanor (ref. 9 )- If 
Treanor js method is applied to the rep- 
resentative equation (7), its stabil- 
ity bounds would be those illustrated 
in sketch (g). The parameter P can 
either be fixed , 2 or varied as the 
integration proceeds. Maximum stabil- 
ity for all crH in the range 
0 > OK > is found when Ph = 8 

and gives johjc = a considerable 
improvement over the Runge-Kutta meth- 
ods shown in sketch (c). For higher 
values of Ph stability can be 
attained for larger -tfH, provided 
they lie in the stability corridor 
shown in the sketch. 



Sketch (g).- Real stability boundaries for 
Treanor* s method. 


HIGHLY STABLE IMPLICIT METHODS 


The explicit methods just discussed have the advantage that they can be 
applied directly to equations (2) and have the theoretical capability of inte- 
grating those equations with arbitrary accuracy. There are many practical 
cases , however, for which their limited stability boundary makes them very 
costly to apply. In such instances the use of implicit methods may be advis- 
able. To compare the efficiency (measured with regard to machine running time) 
of implicit methods with explicit ones is not a simple matter because the 
former require the step by step conversion of equations (2) to equations (3); 
and the subsequent solution of simultaneous algebraic equations. An attempt 
at a rational comparison is presented in reference 8. In this report the capa- 
bilities of the implicit methods are considered without regard to their effi- 
ciency. With these words of caution, we return to the use of the calculation 
step size, h, as the reference parameter in estimates of stability and 
accuracy. 

From the point of view of the computer, the principal difference between 
explicit and implicit methods is that the latter require the solution of M 
simultaneous algebraic equations at each step. From a more abstract point of 
view, the fundamental difference between the two is that the characteristic 
polynomial for the implicit methods is no longer monic, and can always be 
written in the form 


2 Treanor takes P at each step to be the ratio of certain terms in the 
first two steps of a standard fourth-order Runge-Kutta process. For an 
analysis of Treanor *s method see reference 1. 
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P k (orli)E k + P k _ i (cJh)E k ' 1 + • • • + P Q (ah) = 0 


(32) 


where the Pj(oh) are polynomials in Oh. This difference is extremely 
important . 

Consider, for example, the implicit modified Euler method. Let 
represent 

U n+1 = % + \ b(un +1 + Un) 


(33) 


Then Lj_(u* = Ou + f) results in the difference equation 


(l - f on) V: “ (l + | on) % + | h(f n+1 + fn) 


having the characteristic polynomial 


1 - | tfh ) E 


1 + | oh ) = 0 


(3*0 


The single root 


Ai = 


1 + | oh 
1 - | ah 


(35) 


3 - 
2 - 



. 

Unstable 

| 1 m 

4 8 crh 

'Asymptotic to - 

■ 1 
-2 - 

-3 - 

; 


has a magnitude that is less than 1 
(see sketch (h) ) for all Re(tfh) < 0, 
and the method is, as is well known, 
unconditionally stable. 

Not all implicit methods are uncon- 
ditionally stable. For example, if L-j_ 
represents the Adams -Moulton, two-step- 
corrector used implicitly 


u n +2 - u n+i + p2 k(5 u n+2 + 8 u n+i “ u n) 

(36) 

the operation Lp(u* = au + f) results in the characteristic polynomial 


Sketch (h).- Real stability boundary for implicit 
modified Euler method (eq.. ( 33 )) • 


_ _ 5 _ 

12 

which is unstable if ah < -6. 


(l - i ah) E 2 - (l + j| Oh) E + ^ ah - 0 


(37) 
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On the other hand, unconditionally stable methods with any order of 
accuracy can he constructed. The graphic description of 'a class of them is 
shown in sketch (i) for equispaced data. These are sometimes referred to as 

the backward difference formulas and 


<£ 


n- k 


n- k + 1 



Sketch (i).- Unconditionally stable implicit 
methods . 


the first five of them are given in 
reference 10, pages 96-98* These 
methods can be "more stable" than the 
modified Euler method in the sense 
that the roots to their characteristic 
polynomials are smaller for large | crh | . 
For example, the two roots derived 
from the use of 


u n+2 = | ( H1+1 “ %l + 2huA +2 ) 


(38) 


are 


Ai 


2 W 1 + 2qh 
3 - 2ah 


(39a) 


2 - n/ 1 + 2crh 
3 - 2ah 


(39b) 


Both have magnitudes less than 1 for all Re(ah) <0, but behave asymptotically 

as l/J I ahl , rather than as 1 for the modified Euler case. Their variation for 
real ah is shown in sketch (j). Equation (38) has been used in boundary- layer 

studies (ref. 11 ) , and equation (33) 
is extremely popular (ref. 12) for use 
in studying parabolic partial differ- 
ential equations, in which discipline 
they are known as the Crank -Nichols on 
equations . 


Dahlquist has shown (ref. 13 ) 
that only one Lp with p a k + 1 
is unconditionally stable, and in fact, 
the method is that given by equa- 
tion (33)- In this statement p 
refers to the highest order of the 
embedded Taylor series that accurately 
represents e 0 ^ in the expansion of 
Ai = gi(ah) (i.e., p = l - 1 in 
eq. (l4b)), and k is the step number in an equispaced, multistep, linear 
method. Dahlquist is sometimes misquoted as saying that only one implicit 
method is unconditionally stable, namely, equation (33) which has an error 
0(h 3 ). However, his theorem states nothing about the order of accuracy that 
can be achieved in unconditionally stable implicit methods. What it does 
state is that if, for example, an unconditionally stable, four-step, equispaced. 



crh 


Sketch (j).- Heal stability boundary for implicit 
method given by equation (38). 




Ik 



linear , implicit method is used, its embedded Taylor series expansion can be 
accurate, at most, through the term (l/24) (ah) 4 . 

Aside from the hypothesis of equispaced steps in Li, there is another 
qualification underlying Dahlquist T s theorem discussed in the previous 
paragraph; namely, that [A n ] is used only to the first power. Let us consider 
this next in more detail. 

The operator Li converts the (locally) linear differential equations 
into a set of M algebraic equations that must be solved by the methods of 
linear algebra.^ For example, if Li refers to equation (33)^ 

L^(w = [A n ] w + fn) results in 

w n+i = v n + 2 k(Il^nl w n+i + f n + [A n ]v n + f n ) (4o) 

which can be rearranged to form 

(^[l] “2 k[A n ]^ ("Wn+i - w n ) = h[A n ]w n + hf n (4l) 

how it is clear that equation (33) is a special difference -differential 
relationship from the class 

u n+i = u n + b(aiUn + a 2 ui +1 ) + h 2 (b!Un + b 2 Un+i) + • • • (42a) 

which is, in turn, a special case of multistep groups (note, eq. (42a) only 
relates data between one step) with higher order derivatives. In this report, 
however, attention is restricted to equations (42a) and, in particular, to 

u n +i = u n + hau n + ha 2 u n -fi + bsu^+i (42b) 

since it displays all the essential points to be made. Now if u T = ou + f, 
where a is a constant and f is also a constant (i.e., autonomous equations , 
the only kind we will consider in the following ) , u M - o'u f = cr 2 u + erf. Let 
L^ represent equation (42b). Then the operation Lj_(u f = ou + f) results in 
the difference equation 


u n+ 1 = (l + acrh)u I1 + ah(a 2 + b 2 ah)u n + 1 + hf(a + a 2 + ahb 2 ) (43) 

which has the characteristic polynomial 

(l - a 2 crh - b 2 o 2 h 2 )E - (l + anh) = 0 (44) 

The single root expands to give 
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Ax = gi(ah) = ^ _ ~~^ L -°b 2 a a h 2 ~ ^ + ( a + a 2)c7li + [b 2 + a 2 (a + a 2 )]a 2 h 2 

+ [a§(a + a 2 ) + ba(a + 2a 2 )]a 3 h 3 + • • • 


(45) 


Ins isting that Ax represent the expanded e 0 * 1 through the term (l/2)a 2 h 2 
gives 


a 2 = 1 - a 


h 2 = a - | 


(46) 


for which equation (45) reduces to 


Ax = 


1 + ash 


1 + (a - l)cfh +(^| - a)cr 


■ 2 h 2 


= 1 + ah + ~ a 2 h 2 + ^ o 3 h 3 + • • 


(47) 


The exact solution to the difference equation (43) combined with equa- 
tions (46) is 


u n c o | 


1 + aah 


n 


1 + (a - l)ah + (jj - a) O 


2 h 2 


f 

a 


(48a) 


and the exact solution of u' = Ou + f, with constant a and f, is 


- / ah\n f 
u n = c o( e ) - „ 


(48b) 


The particular solution is calculated exactly and the complementary solution 
has an error 0(h 3 ), consistent with the error in equations (3) • 

There remains the discussion of stability for the method Lq defined by 
equations (42b) and (46). Let us first consider the real stability boundary 
in detail and later present the results for complex and imaginary a • When 
a =■ l/3 .> the error is minimized (the method is then 0(h 4 )) and Ax has the 
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Sketch (k) . - Real stability boundary for method 
given by equation ( 52 b) • 


variation shown in sketch (k) . This 
form of the method is not only uncondi- 
tionally stable, it is also stable for 
all positive Oh £ 6. In general, the 
value of Ai(crh) given by equation (47) 
crosses the line +1 when 


ah = 0 

Oh = 2/(1 - 2a) 
and the line -1 when 


( 49 ) 


t 


x, 



( 50 ) 



Sketch (Z).- Real stability boundary for method 
given by equation ( 52 c). 


Thus, between - 3/2 £ a ^ l /2 the value 
of Aq is never less than -1 for any 
Ph, and is above the line 1 only in the 
range 0 < ah < 2/(l - 2 a). Hence, the 
maximum real stability of the method is 
achieved (see sketch (Z)) when 
a = -3/2, in which case the error is 
~( ll/l2) Cah) 3 , and the |Ai| > 1 only 
for 0 < ah <1/2. Of course, the accu- 
racy for ah > l/2 is completely lost, 
but if the positive eigenvalues are 
parasitic -saddle, this is immaterial. 


The complete stability properties 

for all complex ah are shown, for the case a = -3/2, in figure 1 in the 
(A, ah) plane, where the lines ±1 represent the stability boundaries. Figure 2 
shows some of the same results in the complex plane where Re(Ai) is plotted 
against Im(Ai), and the unit circle is the stability boundary. 

When represents 


u n+i " % + htaun + (l - a)u n + 3 _] + ^a - %+! (51) 

the operation Lj_(w = [A n ]w + f n ) results in the difference equation 
([I] + h(a - l)[A n ] + h 2 - a) [A n ] 2 )) (v n+1 - w n ) 

= k (jl] + it (a - [A n ]^ (jA n ]w n + fn) 

(52a) 

IT 



For the most accurate case, a = l/3 and 


^[I] - | h[A n ] + j h p ~[A n ] 2 ^) (v n+1 - w n ) = h ^[1] - j h[A n ]^ ([A n ]w n + f n ) 

( 52 b) 


and for the most stable case, a = - 3/2 and 


[I] - | h[A n ] + 2 h a [A n ] 2 ] (w n+1 - w n ) = h([l] + h[A])([A n ]w n + fn) 


( 52 c) 


Equations (52) are actually programmed when the method is used on -a digital 
computer. They can be said to violate the Dahlquist theorem discussed above; 
that is, equation (42b) is an unconditionally stable method for which (in the 
notation used in presenting the theorem) p = 3 and k = 1. But, as has been 
pointed out, such a violation does not really occur, since Dahlquist's theorem 
is based on the premise that [A n ] is used only to the first power. 


The practical application of methods such as these must, of course, be 
made with care. They are designed for quasi-linear differential equations for 
which the standard general methods are unsatisfactory and something is known, 
or can be determined about, the eigenvalue structure of the associated [A n ] 
matrix. An example of such a case is given by the quasi-linear equation 


(l-w)^j:= 2 w-t + 0.5 


(53) 


with the initial conditions 


to = o 

w 0 = (7 - 5 nT 2)/2 


( 5*0 


It is easy to show that equation (53) has a saddle point (see, e.g., ref. lb) 
at v = 1 and t - 2*5* Introduce the new independent variable such that 


and equation ( 53 ) becomes 



^ = [A]w + f 
ds 


(55) 


( 56 a) 
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where 


Saddle curve 
for which 

W = 1 + (1-2.5) /(I + -/2 - ) vB‘ 


W = 1.0 + C, e 0 " 1 S 


' Saddle point 


^ = [¥,s] 


[A] = 


2 -1 


1 0 


(56b) 


(56c) 


fT = [0.5, 1.0] 


(56d) 


~~ ^8 K6 za 3.2 Notice that the introduction of equa- 

* tion (55) has transformed the nonauton- 

omous equation (53) into the autonomous 
(rn).- Saddle behavior for equation (53)- form (56). The latter may, therefore, 

he integrated using the methods defined 
saddle b y equations (52) . 


Saddle 

asymptote 


The eigenvalues of [A] are 

CTi = - (sj 2 - 1) » -0. 4l4 


Saddle point 


> 


• CTs — + ( v/2* + l) & 2 . 4l4 

I j 1 I l l l 

0 4 8 ' s 2 16 20 24 The initial values are located at 

point (B) in sketch (m) and fall 

Sketch (n).- Locus of compated points for v in exactly on the saddle curve shown. Any 
( w_j s) plane. standard numerical integrating scheme 

applied to equations (56) will diverge 
in typical fashion, following the 

1.2 r curves (BB f ) or (BB n ) shown in the 

Saddle pomU 

sketch. This is Because the coefficient 
of the term with e^sk (which is exactly 

‘ 8 ~ zero along the saddle curve) cannot he 

w e ° * held to zero in finite place arithmetic, 

.4 - . e and an instability will ensue if the 

• 0 " ' numerical method itself is unstable for 

. t ( ( ( ( f positive eigenvalues in this range. 

0* .4 .8 1.2 1.6 2.0 2.4 2.8 

1 Using the method described by equa- 

tion (52c) with a step size h = O.25 

Sketch (o).- Locus of computed points in (t, v) gave thfi resu i ts shown in sketches (n) 

and (o). The use of equation (55) 

transforms the saddle point to infinity in terms of s, the independent vari- 
able used in the numerical integration, so the saddle point is approached 
asymptotically in the (w,s) and (t,s) planes as is shown in sketch (n). (The 
behavior for t(s) was identical.) When w is plotted as a function of x, 
however, the "physical” nature of the variation is regained and shown in 
sketch (o). Notice that the value of crh is about -0.103 for cXih, which 
reduces the local error to about 0.001. The value of the parasitic-saddle 
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eigenvalue is about O. 603 . According to sketch (Z), this gives a value of 
|A(cy 2 h) | < 1 , so the application of the method should be (and was) stable. 
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0 = 1 


Figure 2.- Some results shorn in figure 1 redisplayed in the complex plane 
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